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We present an analytic strong-coupling approach to the phase diagram and elementary excitations 
of the Jaynes-Cummings-Hubbard model describing a superfluid-insulator transition of polaritons 
in an array of coupled QED cavities. In the Mott phase, we find four modes corresponding to parti- 
cle/hole excitations with lower and upper polaritons, respectively. Simple formulas are derived for 
the dispersion and spectral weights within a strong-coupling random-phase approximation (RPA). 
The phase boundary is calculated beyond RPA by including the leading correction due to quantum 
fluctuations. 



The recent experimental success in engineering strong 
interactions between photons and atoms in high-quality 
micro-cavities opens up the possibility to use light- 
matter systems as quantum simulators for many-body 
physics [jy]. A prominent example is the superfluid- 
insulator transition of polaritons in an array of coupled 
QED cavities as described by the Jaynes-Cummings- 
Hubbard model (JCHM) [3, The competition be- 
tween strong atom-photon coupling, giving rise to an ef- 
fective photon repulsion (localization), and the photon 
hopping between cavities (delocalization) leads to a quan- 
tum phase diagram featuring Mott lobes [S] reminiscent 
of those of ultracold atoms in optical lattices as described 
by the seminal Bose-Hubbard model (BHM) [4]. The 
JCHM can be implemented, e^. using single atoms 5], 
excitons ^] , or Cooper pairs [7| . The striking advantage 
of using coupled micro-cavities, with respect to their op- 
tical lattice counterparts, is their individual accessibility. 

Although the coupling between QED cavities has not 
yet been implemented experimentally, the JCHM stimu- 
lated exciting theoretical work over the last three years, 
suggesting that polariton systems can be used to simu- 
late various strongly correlated and exotic phases [1, 0, 
E El [alii. While such exploratory work is justified 
in its own right, still very little is known about the fun- 
damental excitations of the JCHM. The phase boundary 
of the superfluid-insulator transition has been calculated 
using mean-field decoupling 

0,13, Monte-Carlo [H, El 
and variational cluster approaches [3, [13] in two and 
three dimensions. Only two papers have explored the 
fundamental excitations of the system [3, [I3|- 
these results, even on a mean-field level, rely on more 
or less heavy numerical computation due to the intri- 
cate composite nature of polaritons. To the best of our 
knowledge, no analytic results are available neither for 
the phase boundary nor for the excitations of the JCHM. 

In this paper, we show that a linked-cluster expansion 
pioneered for the Fermi-Hubbard model (FHM) [l8| and 
recently applied to the BHM [l^ can be used to obtain 
simple, analytic formulas for the phase diagram as well 
as excitation spectra for arbitrary temperatures, detun- 
ing parameter and lattice geometries. We find two new 
modes, which have been overlooked in previous numerical 



approaches and discuss dispersion and spectral weights 
within strong-coupling RPA. Furthermore, we study the 
effect of quantum fluctuations on the phase boundary of 
the superfluid-insulator transition. In two dimensions, 
our results agree well with recent Monte-Carlo calcula- 
tions. In three dimensions, where numerical data is cur- 
rently not available, we present a quantitatively accurate 
calculation of the quantum phase diagram. 
The Hamiltonian of the JCHM is given by 



H 



a Qi 



(1) 



where hf-^ denotes the local Jaynes-Cummings Hamilto- 
nian hj'-^ — u)ca\ai + L0xafi7~ + g{(J^ai + o-~a|) with 
cavity index i, photon creation (annihilation) operators 
a^^'' and atomic raising (lowering) operators cr^^ \ The 
cavity mode frequency is ujc, the two atomic levels are 
separated by the energy lu^ and the atom-photon cou- 
pling is given by g (we set h ^ 1). The total number 

of excitations, i.e., polaritons TV = J^ii'^l'^i + '^t'^I)^ 
conserved and fixed by the chemical potential /i. The 
third term in ^ describes the delocalization of photons 
over the whole lattice due to hopping between nearest 
neighbour cavities with amplitude J. It competes with 
an effective on-site repulsion between photons mediated 
by the atom-photon coupling. This competition leads to 
Mott lobes in the quantum phase diagram 0] . 

A rough estimate for the size of these Mott lobes can be 
obtained by calculating perturbatively the cost of adding 
or removing a polariton. In the atomic limit (J = 0) 
the eigenstates of the Hamiltonian ([T]) are the dressed 
polariton states labelled by the polariton number n and 
upper/lower branch index cr = ±. For n > they can be 
written as a superposition of a Fock state with n photons 
plus atomic ground state |n, g) and {n — 1) photons with 
the atom in its excited state \{n — 1), e), 

\n+) = sinenlrt,^) -f cos6'„|(7i- l),e) , 

= cos^nlTi,^) - sin^nKn - 1) ,e) , (2) 

with the angle tan6'„ = 2g^/n/{5 + 2xn), Xn = 
\/ g^n + (5^/4 and the detuning parameter 8 — ujc ^ '-Ox- 



The corresponding eigenvalues are 

<^n^ -[t^- ^c)n- 5/2 + (JXn, a = ±. (3) 

The zero polariton state |0— ) = \0,g) is a special case 
with = 0. Upper and lower polariton energies are 
separated by the Rabi splitting ri„ = 2xn- For the cal- 
culation of the quantum phase diagram at small tunnel- 
ing J ^ g we can neglect the upper polariton branch. 
In a straightforward first-order degenerate perturbation 
theory in J, the chemical potentials at which the addi- 
tion/removal of a lower polariton costs no energy is given 
by (we assume a hypercubic lattice in D dimensions) 

Mp-^c - Xn-Xn+i-2DJ{f--,f + Oijyg), 
fih-LUc = Xn-i-Xn + 2DJ{f--f+0{J^/g)-{A) 

the matrix elements f!^^'^ = {na\a^\{n — 1) v) can be eas- 
ily expressed in terms of the angle 0„ and are given by 
/r = (V^ + (Ti^V^r^)/2 for n > 1 (/r = 1/V2) 
at zero detuning ((5 = 0). The two equations in ^ de- 
fine the upper and lower phase boundary in the quan- 
tum phase diagram in Fig. [T] for small values of the hop- 
ping parameter J /g. The point where the two lines meet 
(i.e., jjLp = ^h) is Jc ~ O.lg for n = 1 and represents 
an upper limit for the size of the first Mott lobe. By 
going to higher order in the perturbative expansion for 
the ground-state energy and subsequent resummation of 
the strong-coupling series, one could in principle deter- 
mine the exact location of the phase boundary. This has 
recently been achieved for the BHM (2^, [2l| . 

In this paper, we will follow a different approach and 
study directly the photonic Matsubara Green's func- 
tion Gy (r; t') — — {T ai{T)aj[T')) with the time-ordering 
operator T and the Heisenberg operator aj{T') ~ 
e^"^ a^e"^"^ . A suitable method for the evaluation of 
the Matsubara Green's function is a linked-cluster ex- 
pansion in terms of local cumulants originally devel- 
oped by Metzner et al. [3l for the FHM. Each term 
of the linked-cluster expansion can be written diagram- 
matically in terms of n-particle cumulants represented 
by 2n-leg vertices and tunneling matrix elements sym- 
bolized by propagating lines connecting two vertices. 
The strong-coupling expansion provided by the linked- 
cluster method is applicable to the JCHM because (i) the 
atomic limit Hamiltonian is local, (ii) anomalous averages 
of the photon operator with respect to the eigenstates 
of the local Hamiltonian vanish, i.e., {na\{a')'^\n(T) = 
for k G N (since a single photon excitation always 
changes the polariton number). 

In the atomic limit, the Green's function is given by 
Goij(r;r') = Goi(T; r') (5,^ = -{Tai{T)ai{T'))o Sij, where 
the average (. . . )o is taken with respect to the eigenstates 
of the local Hamiltonian, i.e., the first two terms in ([T]). 
For a spatially homogeneous system we can drop the site 
index i and straightforwardly calculate Gq{t;t'). After 
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FIG. 1: Quantum phase diagram for a hypercubic lattice in 
D — 2 (left figure) and D = 3 (right figure). We show two 
Mott lobes for n — 1,2 at zero detuning 5 = 0. We compare 
first-order perturbation theory (dashed), RPA (dot-dashed) 
and quantum fluctuations (solid) with recent results from a 
quantum Monte-Carlo (filled dots) [ij and variational cluster 
(stars) 17] approach. The two crosses (left figure) connected 
by the dotted line mark the parameter values used in Fig. [5] 
The insets show the critical hopping strength Jc/g at the tip 
of the lobe as a function of the detuning S for n — 1 and n — 2 
calculated within strong-coupling RPA in D = 3 dimensions. 

a Fourier transformation we obtain 

with the partition function Z = X^ncr e^'^'^" and bosonic 
Matsubara frequencies uJm — 27rm//3 (/? = l/iksT) with 
temperature T and Boltzmann constant fee)- In ([5]) we 
defined zi/f = {fl^'^Y and A^'^ = e(^, - e^„i. We sum 
an infinite set of diagrams by calculating the irreducible 
part of the Green's function K{]<i, ujm) which is connected 
to the full Green's function via the equation G(k, Wm) = 
i4r(k, uo,n)l[^ — J(k)ii'(k, Um)] with the lattice dispersion 
J(k) = 2J^^^cosk-ai (a^ denotes a lattice vector). 
To second order in J we get 



= GoM + 2DJ^Q{lo^) (6) 

with = /g^dTdTidT2G(2)(TiO;T2r)Go(T2;ri)e'""^ 

The quantum fluctuation correction Q{uJm) involves 
the two-particle cumulant G*^^-' (tiO; r2T), which is 
related to the local two-particle Green's func- 
tion Gy{TiO;T2T) = (Taj(ri)aj(0)a^(r2)ai(r))o via 
G(2)(TiO;r2T) = G^'^(riO;T2r) - Go(ti; r2)Go(0; r) - 
Go(ti; r)Go(0; r2). The algebraic expressions for the 
two-particle cumulant and the quantum correction 
Q{ujn) are lengthy and will be given elsewhere [2^ . 

The inverse Green's function tells us immediately 
about the phase boundary G^^(0, 0)|j^(^) = and 
the dispersion relation G~^{'k,iujm ^ co + iO~^) — 0. 
We introduce the strong-coupling self-energy E(k, Wm) 
via G(k,w„i)~i = GoiuJm)^^ - S(k,w,„) and ob- 
tain from ^ to second order I](k, cj^) = J(k) -I- 
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2: Particle/hole dispersion of the conventional modes 
(p h) ^'^'^ D = 2 and n = 1 at zero detuning 5 = 0. Left 
figure: Deep inside the Mott insulator (see cross in Fig. [l]) 
for {fj, — iOc)/g — —0.78 and J/g = 0.01. The inset shows the 
conventional modes uj7 , , together with the conversion mode 

Up . Right figure: At the phase boundary (see cross in Fig. [l| 
for (/I — uJc)/g ~ —0.78 and J/g = 0.04. The inset shows 
the bandwidth W of the conventional particle mode ujp as a 
function of detuning S at the tip of the lobe. 

2DJ'^Q{uJm)/Go{uJrn)- The first term on the r.h.s. is 
usually called the strong-coupling random-phase approx- 
imation (RPA) [2^, whereas the second term denotes 
the leading correction due to quantum fluctuations. The 
RPA corresponds to a summation of all self-avoiding 
walks (chain diagrams) through the lattice. The lead- 
ing quantum correction includes in addition all one-time 
forward/backward hopping processes (bubble diagrams) 
between two neighbored sites. Although we have car- 
ried out our calculations at finite temperature, we will 
only consider the T = case from now on. The effect of 
thermal fluctuations will be studied elsewhere [2^ . 

At the quantum phase transition the energy of long 
wavelength fluctuations vanishes and we get an explicit 
expression for the critical hopping strength 



J-i =DGoiO){l + Jl + 2Q{0)/{DGf,m) 



(7) 



If we ignore the second term under the square root in 
d?]) we obtain the RPA phase boundary 1/Jc = 2D Go{0) 
shown as a dashed-dotted line in Fig. [T] We have checked 
that this analytic result agrees exactly with the phase 
boundary obtained from a numerical decoupling mean- 
field approach as presented in 0, Q ■ We can go beyond 
this result by including quantum fluctuations in ([7]) lead- 
ing to an improved phase boundary (solid line in Fig.[T|). 
In two dimensions our result agrees well with recent 
Monte-Carlo calculations 1141 and confirms the smooth- 
ness of the lobes found in [l4| (different from , where a 
variational cluster approach has been used). There is still 
a small deviation near the tip of the lobe, where quan- 
tum corrections are most important. In the right panel of 
Fig. [1] we present quantitatively accurate results for the 
phase diagram in D = 3. Note, that the strong-coupling 
expansion becomes more accurate in higher dimensions. 

We now turn to the discussion of excitations, which we 
calculate within RPA. An analytic continuation of the 
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FIG. 3: Particle/hole gaps (left figure) and effective 

masses (right figure) of the conventional modes for D — 

2 and n = 1 at zero detuning 5 = as a function of the 
tunneling strength J/g for {/j, — ujc)/g = —0.78 (along the 
dotted line in Fig. [IJ . 

Matsubara Green's function via iojn uj + «0+ yields 
the retarded real-time Green's function. Its poles and 
residues yield the dispersion relations and mode strengths 
of the fundamental excitations. In general, we obtain four 
poles, the conventional (Bose-Hubbard like) lower polari- 
ton particle (hole) modes uj'^^ and two modes wj^ 
which correspond to an upper polariton particle (hole) 
excitation (conversion modes). The presence of the lat- 
ter signals a clear deviation from the usual Bose-Hubbard 
like physics and is due to the composite nature of polari- 
tons. The conversion modes exist already in the atomic 
limit and were also found in a recent Monte Carlo study 
[HI , but have been overlooked by a variational cluster ap- 
proach [l7j|. One reason might be that their bandwidth 
and strengths are very small as compared to the conven- 
tional modes. We can thus set their dispersion relations 



^(^pji) and mode strength s^^ 
Lut , , R:! At ,s and st ,^ 

ip,h) {p,h) (p,h) 

limit particle (hole) gaps = A"^^^-^ (A^ = and 

I (z^ = —2^'^), respectively. If 



approximately equal to 
zt i.\ with the atomic- 



mode strengths = 
we neglect their contribution to the one-particle cumu- 
lant (O, we can derive simple analytic formulas for the 
dispersion relations of the conventional modes 



UJ 



(A+ - J(k) z+ ± n) /2 



(8) 



with n = y^A^ 
breviations Aj- 



J(k)2 z\ - 2 J(k) z_ A_ and the ab- 

Ap ± A^ and z± = ± zj^ . The 
strength of the modes are given by 



^+^(p,/i) ^(p,/0^(/i,p) ^(fe,p)^(p,/i) 



(9) 



The dispersions in ([S]) are plotted in Fig. [5] deep inside 
the Mott regime and at the tip of the lobe with n = 1. 
The energy needed for a conventional excitation is an 
order of magnitude smaller than for a conversion excita- 
tion. The strengths of the conventional modes grow 
with increasing tunneling strength (and diverge at the tip 
of the lobe similar to what is found for the BHM ^24] ) 



4 



'(p h) stays approximately constant. If we are only 
interested in low energy excitations we can thus indeed 
neglect contributions from upper polaritons. 

Deep inside the Mott insulator particle and hole exci- 
tations are gapped. If the phase boundary is approached 
away from the tip of the lobe, either the particle (upper 
phase boundary) or the hole (lower phase boundary) gap 
vanishes linearily A ~ | J — Jc|. At the phase boundary 
the dispersion relation remains quadratic u ^ k^. The 
situation changes at the tip of the lobe, where particle 
and hole gaps vanish simultaneously with a square-root 
behavior A ^ \ J — Jc\^^^, while their dispersions become 
linear ~ fc (see Fig. [2]) . This indicates a special transi- 
tion at the tip of the lobe, where also the effective masses 
vanish (see Fig. [3]), reminiscent of an emergent particle- 
hole symmetry. In the opposite limit J ^ the effective 
masses approach infinity as = ±1/(2 J z'^p i^^)- The 

behavior of the dispersion at the critical point J = Jc in 
the long- wavelength limit k — > determines the dynami- 
cal critical exponent z defined hy lu ~ with the 
diverging correlation length ^ ~ | J — Jc| " and its as- 
sociated critical exponent v. From the discussion above, 
we conclude that the dynamical critical exponent has the 
generic value z = 2 everywhere in the phase diagram ex- 
cept for the special multi-critical point at the tip of the 
lobe where it changes to z = 1. Very recently this result 
has been confirmed using an effective action approach 
[2^. At k = the gap vanishes as A ~ | J — Jd^" when 
the tunneling strength approaches its critical value Jc- 
Thus we have v = 1/2 everywhere in the phase diagram. 
We conclude that at least on a mean-field RPA level, the 
JCHM has the same critical exponents as the BHM ji| 
and is thus in the same universality class. The results 
in ([8]) and ([9]) can be seen as a generalization of the ex- 
pressions for the dispersion relations and mode strengths 
of the BHM ^ . The usual BHM physics is retrieved 
if one ignores the upper polariton mode in ([5|) and sets 
fn=fn = V^ and = e„ = {U/2)n{n - 1) - fin. 

The physics of the JCHM is richer due to the pres- 
ence of the experimentally important detuning param- 
eter 5. Indeed, detuning from the resonance (5 = al- 
lows to fine-tune the effective repulsion between photons 
and drive the system from the Mott into the superfluid 
state [2, 14|. In Fig. [1] we show that the critical hopping 
strength decreases as a function of the detuning parame- 
ter S independent of its sign for n = 2. This is generally 
true for Mott lobes with more than one polariton n > 1. 
However, the n = 1 Mott lobe becomes smaller for nega- 
tive and larger for positive detuning. This distinct behav- 
ior is due the special nature of hole excitations from the 
n = 1 polariton state to the zero-polariton state, which 
has an energy eg = independent of detuning. 

At the end we comment on the experimental verifi- 
ability of our results. A setup based on quantum dot 



excitons embedded in a photonic crystal was discussed 
in ^27]. After proper initialization of a Mott insulator 
ground-state, the dispersion relations of the excitations 
could be measured using transmission spectroscopy [28| 
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